// @HEADER
//*********************************************************************//
//  SiYuan: A numerical PDE solver                                     //
//  Copyright (2022) YUAN Xi                                           //
//  This Software is released under the BSD 2-Clause license detailed  //
//  in the file "LICENSE" in the top-level SiYuan directory            //
//*********************************************************************//
// @HEADER

#ifndef _MODEL_EVALUATOR_FACTORY_IMPL_HPP
#define _MODEL_EVALUATOR_FACTORY_IMPL_HPP

#include "Thyra_ModelEvaluator.hpp"
#include "Teuchos_Assert.hpp"
#include "Teuchos_as.hpp"
#include "Teuchos_DefaultMpiComm.hpp"
#include "Teuchos_AbstractFactoryStd.hpp"

// Piro solver objects
#include "Thyra_NonlinearSolver_NOX.hpp"
#include "Piro_ConfigDefs.hpp"
#include "Piro_NOXSolver.hpp"
#include "Piro_LOCASolver.hpp"
#include "Piro_RythmosSolver.hpp"
#include "Piro_TempusSolver.hpp"
#include "Piro_SolverFactory.hpp"
#include "LOCA_Thyra_SaveDataStrategy.H"

#include "PanzerAdaptersSTK_config.hpp"
#include "Panzer_Traits.hpp"
#include "Panzer_GlobalData.hpp"
#include "Panzer_BC.hpp"
#include "Panzer_BasisIRLayout.hpp"
#include "Panzer_DOFManager.hpp"
#include "Panzer_DOFManagerFactory.hpp"
#include "Panzer_BlockedDOFManager.hpp"
#include "Panzer_BlockedDOFManagerFactory.hpp"
#include "Panzer_LinearObjFactory.hpp"
#include "Panzer_TpetraLinearObjFactory.hpp"
#include "Panzer_ThyraObjContainer.hpp"
#include "Panzer_BlockedTpetraLinearObjFactory.hpp"
#include "Panzer_ModelEvaluator.hpp"
#include "Panzer_WorksetContainer.hpp"
#include "Panzer_String_Utilities.hpp"
#include "Panzer_GlobalIndexer_Utilities.hpp"
#include "Panzer_ExplicitModelEvaluator.hpp"
#include "Panzer_ParameterLibraryUtilities.hpp"

#include "Panzer_STK_Interface.hpp"
#include "Panzer_STK_ExodusReaderFactory.hpp"
#include "Panzer_STK_LineMeshFactory.hpp"
#include "Panzer_STK_SquareQuadMeshFactory.hpp"
#include "Panzer_STK_SquareTriMeshFactory.hpp"
#include "Panzer_STK_CubeHexMeshFactory.hpp"
#include "Panzer_STK_CubeTetMeshFactory.hpp"
#include "Panzer_STK_MultiBlockMeshFactory.hpp"
#include "Panzer_STK_CustomMeshFactory.hpp"
#include "Panzer_STK_SetupUtilities.hpp"
#include "TianXin_STK_Utilities.hpp"
#include "Panzer_STK_WorksetFactory.hpp"
#include "Panzer_STKConnManager.hpp"
#include "Panzer_STK_ParameterListCallback.hpp"
#include "Panzer_STK_ParameterListCallbackBlocked.hpp"
#include "Panzer_STK_IOClosureModel_Factory_TemplateBuilder.hpp"
#include "Panzer_STK_ResponseEvaluatorFactory_SolutionWriter.hpp"
#include "Panzer_STK_SetupLOWSFactory.hpp"
#include "Panzer_ResponseEvaluatorFactory_Functional.hpp"

#include "ResponseEvaluatorFactory.hpp"

#include "Tempus_StepperFactory.hpp"
#include "TempusObserverFactory.hpp"
#include "RythmosObserverFactory.hpp"
#include "NOXObserverFactory.hpp"

#include "InitialCondition.hpp"
#include "Eigensolver.hpp"
#include "ROLSolver.hpp"

#include <vector>
#include <iostream>
#include <fstream>
#include <cstdlib> // for std::getenv
#include <random>

namespace SiYuan {

  template<typename ScalarT>
  void ModelEvaluatorFactory<ScalarT>::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
  {
    paramList->validateParametersAndSetDefaults(*this->getValidParameters());
    
    // add in some addtional defaults that are hard to validate externally (this is because of the "disableRecursiveValidation" calls)
  //  if(!paramList->sublist("Initial Conditions").isType<bool>("Zero Initial Conditions"))
  //    paramList->sublist("Initial Conditions").set<bool>("Zero Initial Conditions",false);
  //  paramList->sublist("Initial Conditions").sublist("Vector File").validateParametersAndSetDefaults(
  //    getValidParameters()->sublist("Initial Conditions").sublist("Vector File"));

    this->setMyParamList(paramList);
  }

  template<typename ScalarT>
  Teuchos::RCP<const Teuchos::ParameterList> ModelEvaluatorFactory<ScalarT>::getValidParameters() const
  {
    static Teuchos::RCP<const Teuchos::ParameterList> validPL;
    if (is_null(validPL)) {
      Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::rcp(new Teuchos::ParameterList());

      pl->sublist("Physics Blocks").disableRecursiveValidation();
      pl->sublist("Closure Models").disableRecursiveValidation();
      pl->sublist("Responses").disableRecursiveValidation();
      pl->sublist("fluxResponses").disableRecursiveValidation();
      pl->sublist("Boundary Conditions").disableRecursiveValidation();
      pl->sublist("Dirichlet Conditions").disableRecursiveValidation();
      pl->sublist("Neumann Conditions").disableRecursiveValidation();
      pl->sublist("Concentrated Loads").disableRecursiveValidation();
      pl->sublist("Material").disableRecursiveValidation();
      pl->sublist("Solution Control").disableRecursiveValidation();
      pl->sublist("Linear Solver").disableRecursiveValidation();
      pl->set<bool>("Use Discrete Adjoint",false);

      pl->sublist("Mesh").disableRecursiveValidation();

    /*  pl->sublist("Initial Conditions").set<bool>("Zero Initial Conditions",true);
      pl->sublist("Initial Conditions").sublist("Transient Parameters").disableRecursiveValidation();
      pl->sublist("Initial Conditions").sublist("Vector File");
      pl->sublist("Initial Conditions").sublist("Vector File").set("File Name","");
      pl->sublist("Initial Conditions").sublist("Vector File").set<bool>("Enabled",false);*/
      pl->sublist("Initial Conditions").disableRecursiveValidation();

      pl->sublist("Visual Output").set("File Name","SiYuan.exo");
      pl->sublist("Visual Output").set("Write to Exodus",true);
      pl->sublist("Visual Output").set("Output Frequency",1);
      pl->sublist("Visual Output").sublist("Cell Average Quantities").disableRecursiveValidation();
      pl->sublist("Visual Output").sublist("Cell Quantities").disableRecursiveValidation();
      pl->sublist("Visual Output").sublist("Cell Average Vectors").disableRecursiveValidation();
      pl->sublist("Visual Output").sublist("Nodal Quantities").disableRecursiveValidation();
      pl->sublist("Visual Output").sublist("Allocate Nodal Quantities").disableRecursiveValidation();

      // Assembly sublist
      {
        Teuchos::ParameterList& p = pl->sublist("Assembly");
        p.set<int>("Workset Size", 1);
        p.set<int>("Default Integration Order",-1);
        p.set<std::string>("Field Order","");
        p.set<std::string>("Auxiliary Field Order","");
        p.set<bool>("Load Balance DOFs",false);
        p.set<bool>("Lump Explicit Mass",false);
        p.set<bool>("Constant Mass Matrix",true);
        p.set<bool>("Apply Mass Matrix Inverse in Explicit Evaluator",true);
        p.set<bool>("Use Conservative IMEX",false);
        p.set<bool>("Compute Real Time Derivative",false);
        p.set<bool>("Use Time Derivative in Explicit Model",false);
        p.set<bool>("Compute Time Derivative at Time Step",false);
        p.set<Teuchos::RCP<const panzer::EquationSetFactory> >("Equation Set Factory", Teuchos::null);
        p.set<Teuchos::RCP<const panzer::ClosureModelFactory_TemplateManager<panzer::Traits> > >("Closure Model Factory", Teuchos::null);
        p.set<Teuchos::RCP<const panzer::BCStrategyFactory> >("BC Factory",Teuchos::null);
        p.set<std::string>("Excluded Blocks","");
        p.sublist("ALE").disableRecursiveValidation();
      }

      pl->sublist("Block ID to Physics ID Mapping").disableRecursiveValidation();
      pl->sublist("Block ID to Material ID Mapping").disableRecursiveValidation();
      pl->sublist("Debug Options").disableRecursiveValidation();
      pl->sublist("Solver Factories").disableRecursiveValidation();
      pl->sublist("Active Parameters").disableRecursiveValidation();
      pl->sublist("Controls").disableRecursiveValidation();
      pl->sublist("ALE").disableRecursiveValidation(); // this sucks
      pl->sublist("User Data").disableRecursiveValidation();
      pl->sublist("User Data").sublist("Panzer Data").disableRecursiveValidation();
      pl->sublist("Global Data").disableRecursiveValidation();

      validPL = pl;
    }
    return validPL;
  }

  namespace {
    bool hasInterfaceCondition(const std::vector<panzer::BC>& bcs)
    {
      for (std::vector<panzer::BC>::const_iterator bcit = bcs.begin(); bcit != bcs.end(); ++bcit)
        if (bcit->bcType() == panzer::BCT_Interface)
          return true;
      return false;
    }

    Teuchos::RCP<panzer_stk::STKConnManager>
    getSTKConnManager(const Teuchos::RCP<panzer::ConnManager>& conn_mgr)
    {
      const Teuchos::RCP<panzer_stk::STKConnManager> stk_conn_mgr =
        Teuchos::rcp_dynamic_cast<panzer_stk::STKConnManager>(conn_mgr);
      TEUCHOS_TEST_FOR_EXCEPTION(stk_conn_mgr.is_null(), std::logic_error,
                                 "There are interface conditions, but the connection manager"
                                 " does not support the necessary connections.");
      return stk_conn_mgr;
    }

    void buildInterfaceConnections(const std::vector<panzer::BC>& bcs,
                                   const Teuchos::RCP<panzer::ConnManager>& conn_mgr)
    {
      const Teuchos::RCP<panzer_stk::STKConnManager> stk_conn_mgr = getSTKConnManager(conn_mgr);
      for (std::vector<panzer::BC>::const_iterator bcit = bcs.begin(); bcit != bcs.end(); ++bcit)
        if (bcit->bcType() == panzer::BCT_Interface)
          stk_conn_mgr->associateElementsInSideset(bcit->sidesetID());
    }

    void checkInterfaceConnections(const Teuchos::RCP<panzer::ConnManager>& conn_mgr,
                                   const Teuchos::RCP<Teuchos::Comm<int> >& comm)
    {
      const Teuchos::RCP<panzer_stk::STKConnManager> stk_conn_mgr = getSTKConnManager(conn_mgr);
      std::vector<std::string> sidesets = stk_conn_mgr->checkAssociateElementsInSidesets(*comm);
      if ( ! sidesets.empty()) {
        std::stringstream ss;
        ss << "Sideset IDs";
        for (std::size_t i = 0; i < sidesets.size(); ++i)
          ss << " " << sidesets[i];
        ss << " did not yield associations, but these sidesets correspond to BCT_Interface BCs.";
        TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, ss.str());
      }
    }
  } // namespace

  template<typename ScalarT>
  void  ModelEvaluatorFactory<ScalarT>::buildObjects(const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
                                                     const Teuchos::RCP<panzer::GlobalData>& global_data,
                                                     const Teuchos::RCP<const panzer::EquationSetFactory>& eqset_factory,
                                                     const panzer::BCStrategyFactory & bc_factory,
                                                     const panzer::ClosureModelFactory_TemplateManager<panzer::Traits> & user_cm_factory,
                                                     bool meConstructionOn)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(this->getParameterList()), std::runtime_error,
                       "ParameterList must be set before objects can be built!");

    TEUCHOS_ASSERT(nonnull(comm));
    TEUCHOS_ASSERT(nonnull(global_data));
    TEUCHOS_ASSERT(nonnull(global_data->os));
    TEUCHOS_ASSERT(nonnull(global_data->pl));

    // begin at the beginning...
    m_global_data = global_data;

    ////////////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////////////////
    // Parse input file, setup parameters
    ////////////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////////////////

    // this function will need to be broken up eventually and probably
    // have parts moved back into panzer.  Just need to get something
    // running.

    Teuchos::ParameterList& p = *this->getNonconstParameterList();

    // "parse" parameter list
    Teuchos::ParameterList & mesh_params     = p.sublist("Mesh");
    Teuchos::ParameterList & assembly_params = p.sublist("Assembly");
    Teuchos::ParameterList & material_params = p.sublist("Material");
    Teuchos::ParameterList & solncntl_params = p.sublist("Solution Control");
    Teuchos::ParameterList & closure_params = p.sublist("Closure Models");
    Teuchos::ParameterList & output_list = p.sublist("Visual Output");
    Teuchos::ParameterList & dirichlet_list = p.sublist("Dirichlet Conditions");
    Teuchos::ParameterList & neumann_list = p.sublist("Neumann Conditions");
    Teuchos::ParameterList & cload_list = p.sublist("Concentrated Loads");

    Teuchos::ParameterList & user_data_params = p.sublist("User Data");
    Teuchos::ParameterList & panzer_data_params = user_data_params.sublist("Panzer Data");

    Teuchos::RCP<Teuchos::ParameterList> physics_block_plist = Teuchos::sublist(this->getMyNonconstParamList(),"Physics Blocks");

    // extract assembly information
    std::size_t workset_size = Teuchos::as<std::size_t>(assembly_params.get<int>("Workset Size"));
    std::string field_order  = assembly_params.get<std::string>("Field Order"); // control nodal ordering of unknown
                                                                                   // global IDs in linear system
    bool use_load_balance = assembly_params.get<bool>("Load Balance DOFs");

    // this is weird...we are accessing the solution control to determine if things are transient
    // it is backwards!
    m_is_transient  = false;
    m_is_dotdot = false;
    // for pseudo-transient, we need to enable transient solver support to get time derivatives into fill
    std::string solution = solncntl_params.get<std::string>("Piro Solver");
    if( solution=="Tempus" || solution=="Rythmos" || solution=="Eigen") {
      m_is_transient = true;
    }
    else if ( solution == "NOX") {
      if (solncntl_params.sublist("NOX").get<std::string>("Nonlinear Solver") == "Pseudo-Transient")
        m_is_transient = true;
    }
    // for eigenvalues, we need to enable transient solver support to
    // get time derivatives into generalized eigenvale problem
    else if ( solution == "LOCA") {
      if (solncntl_params.sublist("LOCA").sublist("Stepper").get<bool>("Compute Eigenvalues"))
        m_is_transient = true;
    }

    useDiscreteAdjoint = p.get<bool>("Use Discrete Adjoint");

    ////////////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////////////////
    // Do stuff
    ////////////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////////////////

    Teuchos::FancyOStream& fout = *global_data->os;

    // for convience cast to an MPI comm
    const Teuchos::RCP<const Teuchos::MpiComm<int> > mpi_comm =
      Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);

    // Build mesh factory and uncommitted mesh
    ////////////////////////////////////////////////////////////////////////////////////////

    Teuchos::RCP<panzer_stk::STK_MeshFactory> mesh_factory = this->buildSTKMeshFactory(mesh_params);
    Teuchos::RCP<panzer_stk::STK_Interface> mesh = mesh_factory->buildUncommitedMesh(*(mpi_comm->getRawMpiComm()));
    m_mesh = mesh;
    std::cout << "End of mesh construction!\n";
    m_eqset_factory = eqset_factory;

    // setup the physcs blocks
    ////////////////////////////////////////////////////////////////////////////////////////

    std::vector<Teuchos::RCP<panzer::PhysicsBlock> > physicsBlocks;
    {
      // cell ( block id -> physics id ) mapping
      std::map<std::string,std::string> block_ids_to_physics_ids;
      auto bp = p.sublist("Block ID to Physics ID Mapping", true);
      for (Teuchos::ParameterList::ConstIterator entry = bp.begin(); entry != bp.end(); ++entry) {
        std::string dummy_type;
        block_ids_to_physics_ids[entry->first] = entry->second.getValue(&dummy_type);
      }
      // cell ( block id -> material id ) mapping
      std::map<std::string,std::string> block_ids_to_matl_ids;
      auto bm = p.sublist("Block ID to Material ID Mapping", true);
      for (Teuchos::ParameterList::ConstIterator entry = bm.begin(); entry != bm.end(); ++entry) {
        std::string dummy_type;
        block_ids_to_matl_ids[entry->first] = entry->second.getValue(&dummy_type);
      }
      TEUCHOS_TEST_FOR_EXCEPTION(block_ids_to_matl_ids.empty(), std::runtime_error,
                       "Material not defined!");

      for(auto itr: block_ids_to_physics_ids ) {
        std::string element_block_id = itr.first;
        std::string physics_block_id = itr.second;
        std::string material_id = block_ids_to_matl_ids.at(element_block_id);
        const panzer::CellData volume_cell_data(workset_size, mesh->getCellTopology(element_block_id));     

        // find physics block parameter sublist
        TEUCHOS_TEST_FOR_EXCEPTION(!physics_block_plist->isSublist(physics_block_id),
                            std::runtime_error,
                            "Failed to find physics id: \""
                            << physics_block_id
                            << "\" requested by element block: \""
                            << element_block_id << "\"!");
        Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::sublist(physics_block_plist,physics_block_id);
        pl->setName(physics_block_id);
      
        // the physics block nows how to build and register evaluator with the field manager
        Teuchos::RCP<PhysicsBlock> pb = Teuchos::rcp(
          new PhysicsBlock(pl,
				       element_block_id, 
               material_id,
				       assembly_params.get<int>("Default Integration Order"),
				       volume_cell_data,
				       eqset_factory,
				       global_data,
				       m_is_transient));
        m_is_dotdot = m_is_dotdot || pb->isDotDot();

        // we can have more than one physics block, one per element block
        physicsBlocks.push_back(pb);
      }
      m_physics_blocks = physicsBlocks; // hold onto physics blocks for safe keeping
    }

    // add fields automatically written through the closure model
    ////////////////////////////////////////////////////////////////////////////////////////
    addUserFieldsToMesh(*mesh,output_list);

    // finish building mesh, set required field variables and mesh bulk data
    ////////////////////////////////////////////////////////////////////////////////////////

    try {
       // this throws some exceptions, catch them as neccessary
       this->finalizeMeshConstruction(*mesh_factory,physicsBlocks,*mpi_comm,*mesh);
    } catch(const std::logic_error &ebexp) {
       fout << "*****************************************\n\n";
       fout << "Element block exception, could not finalize the mesh, printing block and sideset information:\n";
       fout.pushTab(3);
       mesh->printMetaData(fout);
       fout.popTab();
       fout << std::endl;

       throw ebexp;
    } catch(const std::logic_error & ssexp) {
       fout << "*****************************************\n\n";
       fout << "Sideset exception, could not finalize the mesh, printing block and sideset information:\n";
       fout.pushTab(3);
       mesh->printMetaData(fout);
       fout.popTab();
       fout << std::endl;

       throw ssexp;
    }

  //  mesh->print(fout);
    if(p.sublist("Visual Output").get<bool>("Write to Exodus"))
      mesh->setupExodusFile(p.sublist("Visual Output").get<std::string>("File Name"));

    // build a workset factory that depends on STK
    ////////////////////////////////////////////////////////////////////////////////////////
    Teuchos::RCP<panzer_stk::WorksetFactory> wkstFactory;
    if(m_user_wkst_factory==Teuchos::null)
       wkstFactory = Teuchos::rcp(new panzer_stk::WorksetFactory()); // build STK workset factory
    else
       wkstFactory = m_user_wkst_factory;

     // set workset factory mesh
     wkstFactory->setMesh(mesh);

    // handle boundary and interface conditions
    ////////////////////////////////////////////////////////////////////////////////////////
    std::vector<panzer::BC> bcs;
    panzer::buildBCs(bcs, p.sublist("Boundary Conditions"), global_data);

    // build the connection manager
    ////////////////////////////////////////////////////////////////////////////////////////
    Teuchos::RCP<panzer::ConnManager> conn_manager = Teuchos::rcp(new panzer_stk::STKConnManager(mesh));
    m_conn_manager = conn_manager;

    // build DOF Manager
    ////////////////////////////////////////////////////////////////////////////////////////

    Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> > linObjFactory;
    Teuchos::RCP<panzer::GlobalIndexer> globalIndexer;

    std::string loadBalanceString = ""; // what is the load balancing information
    bool blockedAssembly = false;

    if( panzer::BlockedDOFManagerFactory::requiresBlocking(field_order) ) {

       // use a blocked DOF manager
       blockedAssembly = true;

       panzer::BlockedDOFManagerFactory globalIndexerFactory;

       Teuchos::RCP<panzer::GlobalIndexer> dofManager
         = globalIndexerFactory.buildGlobalIndexer(mpi_comm->getRawMpiComm(),physicsBlocks,conn_manager,field_order);
       globalIndexer = dofManager;

       Teuchos::RCP<panzer::BlockedTpetraLinearObjFactory<panzer::Traits,double,int,panzer::GlobalOrdinal> > bloLinObjFactory
        = Teuchos::rcp(new panzer::BlockedTpetraLinearObjFactory<panzer::Traits,double,int,panzer::GlobalOrdinal>(mpi_comm,
                                                          Teuchos::rcp_dynamic_cast<panzer::BlockedDOFManager>(dofManager)));

       // parse any explicitly excluded pairs or blocks
       const std::string excludedBlocks = assembly_params.get<std::string>("Excluded Blocks");
       std::vector<std::string> stringPairs;
       panzer::StringTokenizer(stringPairs,excludedBlocks,";",true);
       for(std::size_t i=0;i<stringPairs.size();i++) {
          std::vector<std::string> sPair;
          std::vector<int> iPair;
          panzer::StringTokenizer(sPair,stringPairs[i],",",true);
          panzer::TokensToInts(iPair,sPair);

          TEUCHOS_TEST_FOR_EXCEPTION(iPair.size()!=2,std::logic_error,
                        "Input Error: The correct format for \"Excluded Blocks\" parameter in \"Assembly\" sub list is:\n"
                        "   <int>,<int>; <int>,<int>; ...; <int>,<int>\n"
                        "Failure on string pair " << stringPairs[i] << "!");

          bloLinObjFactory->addExcludedPair(iPair[0],iPair[1]);
       }

       linObjFactory = bloLinObjFactory;

       // build load balancing string for informative output
       loadBalanceString = printUGILoadBalancingInformation(*dofManager);
    }
    else {
       // use a flat DOF manager
       panzer::DOFManagerFactory globalIndexerFactory;
       globalIndexerFactory.setUseTieBreak(use_load_balance);
       Teuchos::RCP<panzer::GlobalIndexer> dofManager
         = globalIndexerFactory.buildGlobalIndexer(mpi_comm->getRawMpiComm(),physicsBlocks,conn_manager,field_order);
       globalIndexer = dofManager;

       TEUCHOS_ASSERT(!useDiscreteAdjoint); // safety check
       linObjFactory = Teuchos::rcp(new panzer::TpetraLinearObjFactory<panzer::Traits,double,int,panzer::GlobalOrdinal>(mpi_comm,dofManager));

       // build load balancing string for informative output
       loadBalanceString = printUGILoadBalancingInformation(*dofManager);
    }
    
    TEUCHOS_ASSERT(globalIndexer!=Teuchos::null);
    TEUCHOS_ASSERT(linObjFactory!=Teuchos::null);
    m_global_indexer = globalIndexer;
    m_lin_obj_factory = linObjFactory;
    m_blockedAssembly = blockedAssembly;

    std::vector< Dirichlet<panzer::Traits::RealType> > dirichlets;
    readDirichletCondition( dirichlets, dirichlet_list, mesh, globalIndexer);
    std::vector< NeumannBoundary > neumanns;
  //  readNeumannBoundary( neumanns, neumann_list, mesh );
    std::vector< CLoad > cloads;
    readCLoad( cloads, cload_list, mesh );

    // print out load balancing information
  //  fout << "Degree of freedom load balancing: " << loadBalanceString << std::endl;

    // build worksets
    //////////////////////////////////////////////////////////////
   
    // build up needs array for workset container
    std::map<std::string,panzer::WorksetNeeds> needs;  
    for(std::size_t i=0;i<physicsBlocks.size();i++)
      needs[physicsBlocks[i]->elementBlockID()] = physicsBlocks[i]->getWorksetNeeds();

    Teuchos::RCP<panzer::WorksetContainer> wkstContainer     // attach it to a workset container (uses lazy evaluation)
       = Teuchos::rcp(new panzer::WorksetContainer(wkstFactory,needs));

    wkstContainer->setWorksetSize(workset_size);
    wkstContainer->setGlobalIndexer(globalIndexer); // set the global indexer so the orientations are evaluated

    m_wkstContainer = wkstContainer;

    // find max number of worksets
    std::size_t max_wksets = 0;
    for(std::size_t pb=0;pb<physicsBlocks.size();pb++) {
      const panzer::WorksetDescriptor wd = panzer::blockDescriptor(physicsBlocks[pb]->elementBlockID());
      Teuchos::RCP< std::vector<panzer::Workset> >works = wkstContainer->generateWorksets(wd);
      max_wksets = std::max(max_wksets,works->size());
    }
    user_data_params.set<std::size_t>("Max Worksets",max_wksets);
    wkstContainer->clear(); 

    // Setup lagrangian type coordinates
    /////////////////////////////////////////////////////////////

    // see if field coordinates are required, if so reset the workset container
    // and set the coordinates to be associated with a field in the mesh
    useDynamicCoordinates_ = false;
    for(std::size_t pb=0;pb<physicsBlocks.size();pb++) {
      if(physicsBlocks[pb]->getCoordinateDOFs().size()>0) {
         mesh->setUseFieldCoordinates(true);
         useDynamicCoordinates_ = true;
         wkstContainer->clear(); // this serves to refresh the worksets
                                 // and put in new coordinates
         break;
      }
    }

    // Add mesh objects to user data to make available to user ctors
    /////////////////////////////////////////////////////////////

    panzer_data_params.set("STK Mesh", mesh);
    panzer_data_params.set("DOF Manager", globalIndexer);
    panzer_data_params.set("Linear Object Factory", linObjFactory);

    // If user requested it, short circuit model construction
    ////////////////////////////////////////////////////////////////////////////////////////

    if(!meConstructionOn) return;

    // Setup active parameters
    /////////////////////////////////////////////////////////////

    std::vector<Teuchos::RCP<Teuchos::Array<std::string> > > p_names;
    std::vector<Teuchos::RCP<Teuchos::Array<double> > > p_values;
    if (p.isSublist("Active Parameters")) {
      Teuchos::ParameterList& active_params = p.sublist("Active Parameters");

      int num_param_vecs = active_params.get<int>("Number of Parameter Vectors",0);
      p_names.resize(num_param_vecs);
      p_values.resize(num_param_vecs);
      for (int i=0; i<num_param_vecs; i++) {
        std::stringstream ss;
        ss << "Parameter Vector " << i;
        Teuchos::ParameterList& pList = active_params.sublist(ss.str());
        int numParameters = pList.get<int>("Number");
        TEUCHOS_TEST_FOR_EXCEPTION(numParameters == 0,
                                   Teuchos::Exceptions::InvalidParameter,
                                   std::endl << "Error!  panzer::ModelEvaluator::ModelEvaluator():  " <<
                                   "Parameter vector " << i << " has zero parameters!" << std::endl);
        p_names[i] =
          Teuchos::rcp(new Teuchos::Array<std::string>(numParameters));
        p_values[i] =
          Teuchos::rcp(new Teuchos::Array<double>(numParameters));
        for (int j=0; j<numParameters; j++) {
          std::stringstream ss2;
          ss2 << "Parameter " << j;
          (*p_names[i])[j] = pList.get<std::string>(ss2.str());
          ss2.str("");

          ss2 << "Initial Value " << j;
          (*p_values[i])[j] = pList.get<double>(ss2.str());

          // this is a band-aid/hack to make sure parameters are registered before they are accessed
          panzer::registerScalarParameter((*p_names[i])[j],*global_data->pl,(*p_values[i])[j]);
        }
      }
    }

    // setup the closure model for automatic writing (during residual/jacobian update)
    ////////////////////////////////////////////////////////////////////////////////////////

    panzer_stk::IOClosureModelFactory_TemplateBuilder<panzer::Traits> io_cm_builder(user_cm_factory,mesh,output_list);
    panzer::ClosureModelFactory_TemplateManager<panzer::Traits> cm_factory;
    cm_factory.buildObjects(io_cm_builder);

    // Options about debug output
    /////////////////////////////////////////////////////////////

    bool write_dot_files = p.sublist("Debug Options").get("Write Volume Assembly Graphs",false);
    std::string dot_file_prefix = p.sublist("Debug Options").get("Volume Assembly Graph Prefix","AssemblyGraph_");
    bool write_fm_files = p.sublist("Debug Options").get("Write Field Manager Files",false);
    std::string fm_file_prefix = p.sublist("Debug Options").get("Field Manager File Prefix","AssemblyGraph_");

    // Allow users to override inputs via runtime env
    {
        auto check_write_dag = std::getenv("PANZER_WRITE_DAG");
        if (check_write_dag != nullptr) {
          write_dot_files = true;
          write_fm_files = true;
        }
    }
    
    // setup field manager build
    /////////////////////////////////////////////////////////////

    Teuchos::RCP<FieldManagerBuilder> fmb;
    {
      fmb = buildFieldManagerBuilder(wkstContainer,physicsBlocks,bcs,dirichlets,neumanns,cloads,*eqset_factory,bc_factory,cm_factory,
                                     user_cm_factory,closure_params,*linObjFactory,material_params,
                                     write_dot_files,dot_file_prefix, write_fm_files,fm_file_prefix);
      if( solncntl_params.get<std::string>("Piro Solver")=="Eigen" )
      {
        std::vector<bool> aet = {false, true, false};   // only Jacobain true
        fmb->setActiveEvaluationTypes( aet );
      }
    }

    // build response library
    /////////////////////////////////////////////////////////////

    m_response_library = Teuchos::rcp(new panzer::ResponseLibrary<panzer::Traits>(wkstContainer,globalIndexer,linObjFactory));

    // Setup solver factory
    /////////////////////////////////////////////////////////////

    Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
          buildLOWSFactory(blockedAssembly,globalIndexer,conn_manager,mesh,mpi_comm);

    // Setup physics model evaluator
    /////////////////////////////////////////////////////////////

    double t_init = 0.0;   // why initial time here ?
    Teuchos::RCP<Thyra::ModelEvaluatorDefaultBase<double> > thyra_me
        = Teuchos::rcp(new panzer::ModelEvaluator<double>
                  (fmb, m_response_library,
                                     linObjFactory,
                                     p_names,
                                     p_values,
                                     lowsFactory,
                                     global_data,
                                     m_is_transient, m_is_dotdot,
                                     t_init) );

    // Setup initial conditions
    /////////////////////////////////////////////////////////////
    if ( p.sublist("Initial Conditions").numParams()>0 ) {

      Teuchos::ParameterList pl = p.sublist("Initial Conditions");
      setInitialConditions(pl, *thyra_me, globalIndexer, linObjFactory, physicsBlocks[0]);

      // Write the IC vector into the STK mesh: use response library
      //////////////////////////////////////////////////////////////////////////
      writeInitialConditions(*thyra_me,physicsBlocks,wkstContainer,globalIndexer,linObjFactory,mesh,user_cm_factory,
                           closure_params,
                           p.sublist("User Data"),workset_size);
    }

    m_stkIOResponseLibrary = Teuchos::rcp(new panzer::ResponseLibrary<panzer::Traits>());

    m_physics_me = thyra_me;
  }

  template<typename ScalarT>
  void ModelEvaluatorFactory<ScalarT>::setInitialConditions(Teuchos::ParameterList& pls, 
          Thyra::ModelEvaluatorDefaultBase<ScalarT> & model,
          Teuchos::RCP<panzer::GlobalIndexer>& globalIndexer,
	        Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> >& lof,
          Teuchos::RCP<panzer::PhysicsBlock>& physicsBlock)
  {
	  Teuchos::RCP<panzer::LinearObjContainer> loc = lof->buildGhostedLinearObjContainer();

    Thyra::ModelEvaluatorBase::InArgs<double> nomValues = model.getNominalValues();
    Teuchos::RCP<Thyra::VectorBase<double> > x_vec = Teuchos::rcp_const_cast<Thyra::VectorBase<double> >(nomValues.get_x());
    Teuchos::RCP<Thyra::VectorBase<double> > x_dot;
    if( m_is_transient ) {
      x_dot = Teuchos::rcp_const_cast<Thyra::VectorBase<double> >(nomValues.get_x_dot());
    }

    if( pls.isParameter("Zero") ) {
      if( pls.get<bool>("Zero") )  {
        Thyra::assign(x_vec.ptr(),0.0);
        if( m_is_transient ) Thyra::assign(x_dot.ptr(),0.0);
      } 
    }
    else if( pls.isParameter("Value") ) {
      auto v = pls.get<double>("Value");
      Thyra::assign(x_vec.ptr(),v);
    }
    else if( pls.isParameter("File Name") ) {
      const std::string filename = pls.get<std::string>("FileName");
      lof->readVector(filename,*loc,panzer::LinearObjContainer::X);
    }
    else {
      Thyra::assign(x_vec.ptr(),0.0);
      if( m_is_transient ) Thyra::assign(x_dot.ptr(),0.0);
      for (Teuchos::ParameterList::ConstIterator pl=pls.begin(); pl != pls.end(); ++pl) {
        TEUCHOS_TEST_FOR_EXCEPTION( !(pl->second.isList()), std::logic_error,
				  "Error - All objects in the Initial Conditions sublist must be sublists!" );
        Teuchos::ParameterList& subpls = Teuchos::getValue<Teuchos::ParameterList>(pl->second);

        std::string fdnamepl = subpls.name();
        auto pos = fdnamepl.rfind("->");
        auto fdname = fdnamepl.substr(pos+2);
        int fnum = globalIndexer->getFieldNum(fdname);
        TEUCHOS_TEST_FOR_EXCEPTION( (fnum==-1), std::logic_error,
				      "Error - Field name not found in Initial Conditions definition!" );

        for (Teuchos::ParameterList::ConstIterator spl=subpls.begin(); spl!= subpls.end(); ++spl) {
          TEUCHOS_TEST_FOR_EXCEPTION( !(spl->second.isList()), std::logic_error,
				      "Error - All objects in the Initial Conditions sublist must be sublists!" );
          Teuchos::ParameterList& fd = Teuchos::getValue<Teuchos::ParameterList>(spl->second);
          std::string ebname = fd.name();
          pos = ebname.rfind("->");
          auto ss = ebname.substr(pos+2);

          std::set<panzer::LocalOrdinal> ldofs;
          if( ss=="Element All" ) {
            globalIndexer->getFieldIndex(fdname,ldofs);
          }
          else {
            globalIndexer->getFieldIndex_ElementBlock(fdname, ss, ldofs);
          }

          std::string method = fd.get<std::string>("Method");
          if( method=="Randomize" ) {
                Teuchos::Array<double> rg = fd.get< Teuchos::Array<double> >("Range");
                std::random_device seed_gen;
                std::default_random_engine engine(seed_gen());
                std::uniform_real_distribution<double> rnd(rg[0], rg[1]);
                for( auto ldof: ldofs ) {
                  double d = rnd(engine);
                  Thyra::set_ele(ldof, d, x_vec.ptr());
                }
          }
          else if( method=="Constant" ) {
                double v = fd.get< double >("Value");
                for( auto ldof: ldofs ) {
                  Thyra::set_ele(ldof, v, x_vec.ptr());
                }
          }
        }
      }
    }
    Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<double> >(loc)->set_x_th(x_vec);
    if( m_is_transient ) Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<double> >(loc)->set_dxdt_th(x_dot);
  }

  template<typename ScalarT>
  void ModelEvaluatorFactory<ScalarT>::
  addUserFieldsToMesh(panzer_stk::STK_Interface & mesh,const Teuchos::ParameterList & output_list) const
  {
    // register cell averaged scalar fields
    const Teuchos::ParameterList & cellAvgQuants = output_list.sublist("Cell Average Quantities");
    for(Teuchos::ParameterList::ConstIterator itr=cellAvgQuants.begin();
        itr!=cellAvgQuants.end();++itr) {
       const std::string & blockId = itr->first;
       const std::string & fields = Teuchos::any_cast<std::string>(itr->second.getAny());
       std::vector<std::string> tokens;

       // break up comma seperated fields
       panzer::StringTokenizer(tokens,fields,",",true);

       for(std::size_t i=0;i<tokens.size();i++)
          mesh.addCellField(tokens[i],blockId);
    }

    // register cell averaged components of vector fields
    // just allocate space for the fields here. The actual calculation and writing
    // are done by panzer_stk::ScatterCellAvgVector.
    const Teuchos::ParameterList & cellAvgVectors = output_list.sublist("Cell Average Vectors");
    for(Teuchos::ParameterList::ConstIterator itr = cellAvgVectors.begin();
        itr != cellAvgVectors.end(); ++itr) {
       const std::string & blockId = itr->first;
       const std::string & fields = Teuchos::any_cast<std::string>(itr->second.getAny());
       std::vector<std::string> tokens;

       // break up comma seperated fields
       panzer::StringTokenizer(tokens,fields,",",true);

       for(std::size_t i = 0; i < tokens.size(); i++) {
          std::string d_mod[3] = {"X","Y","Z"};
          for(std::size_t d = 0; d < mesh.getDimension(); d++)
              mesh.addCellField(tokens[i]+d_mod[d],blockId);
       }
    }

    // register cell quantities
    const Teuchos::ParameterList & cellQuants = output_list.sublist("Cell Quantities");
    for(Teuchos::ParameterList::ConstIterator itr=cellQuants.begin();
        itr!=cellQuants.end();++itr) {
       const std::string & blockId = itr->first;
       const std::string & fields = Teuchos::any_cast<std::string>(itr->second.getAny());
       std::vector<std::string> tokens;

       // break up comma seperated fields
       panzer::StringTokenizer(tokens,fields,",",true);

       for(std::size_t i=0;i<tokens.size();i++)
          mesh.addCellField(tokens[i],blockId);
    }

    // register ndoal quantities
    const Teuchos::ParameterList & nodalQuants = output_list.sublist("Nodal Quantities");
    for(Teuchos::ParameterList::ConstIterator itr=nodalQuants.begin();
        itr!=nodalQuants.end();++itr) {
       const std::string & blockId = itr->first;
       const std::string & fields = Teuchos::any_cast<std::string>(itr->second.getAny());
       std::vector<std::string> tokens;

       // break up comma seperated fields
       panzer::StringTokenizer(tokens,fields,",",true);

       for(std::size_t i=0;i<tokens.size();i++)
          mesh.addSolutionField(tokens[i],blockId);
    }

    const Teuchos::ParameterList & allocNodalQuants = output_list.sublist("Allocate Nodal Quantities");
    for(Teuchos::ParameterList::ConstIterator itr=allocNodalQuants.begin();
        itr!=allocNodalQuants.end();++itr) {
       const std::string & blockId = itr->first;
       const std::string & fields = Teuchos::any_cast<std::string>(itr->second.getAny());
       std::vector<std::string> tokens;

       // break up comma seperated fields
       panzer::StringTokenizer(tokens,fields,",",true);

       for(std::size_t i=0;i<tokens.size();i++)
          mesh.addSolutionField(tokens[i],blockId);
    }
  }

  template<typename ScalarT>
  void ModelEvaluatorFactory<ScalarT>::
  writeInitialConditions(const Thyra::ModelEvaluator<ScalarT> & model,
                         const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
                         const Teuchos::RCP<panzer::WorksetContainer> & wc,
                         const Teuchos::RCP<const panzer::GlobalIndexer> & ugi,
                         const Teuchos::RCP<const panzer::LinearObjFactory<panzer::Traits> > & lof,
                         const Teuchos::RCP<panzer_stk::STK_Interface> & mesh,
                         const panzer::ClosureModelFactory_TemplateManager<panzer::Traits> & cm_factory,
                         const Teuchos::ParameterList & closure_model_pl,
                         const Teuchos::ParameterList & user_data_pl,
                         int workset_size) const
  {
    Teuchos::RCP<panzer::LinearObjContainer> loc = lof->buildLinearObjContainer();
    Teuchos::RCP<panzer::ThyraObjContainer<double> > tloc = Teuchos::rcp_dynamic_cast<panzer::ThyraObjContainer<double> >(loc);
    tloc->set_x_th(Teuchos::rcp_const_cast<Thyra::VectorBase<double> >(model.getNominalValues().get_x()));

    Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > solnWriter
        = initializeSolnWriterResponseLibrary(wc,ugi,lof,mesh);

    {
       Teuchos::ParameterList user_data(user_data_pl);
       user_data.set<int>("Workset Size",workset_size);

       finalizeSolnWriterResponseLibrary(*solnWriter,physicsBlocks,cm_factory,closure_model_pl,workset_size,user_data);
    }

    // initialize the assembly container
    panzer::AssemblyEngineInArgs ae_inargs;
    ae_inargs.container_ = loc;
    ae_inargs.ghostedContainer_ = lof->buildGhostedLinearObjContainer();
    ae_inargs.alpha = 0.0;
    ae_inargs.beta = 1.0;
    ae_inargs.evaluate_transient_terms = false;

    // initialize the ghosted container
    lof->initializeGhostedContainer(panzer::LinearObjContainer::X,*ae_inargs.ghostedContainer_);

    // do import
    lof->globalToGhostContainer(*ae_inargs.container_,*ae_inargs.ghostedContainer_,panzer::LinearObjContainer::X);

    // fill STK mesh objects
    solnWriter->addResponsesToInArgs<panzer::Traits::Residual>(ae_inargs);
    solnWriter->evaluate<panzer::Traits::Residual>(ae_inargs);
  }

  //! build STK mesh from a mesh parameter list
  template<typename ScalarT>
  Teuchos::RCP<panzer_stk::STK_MeshFactory> ModelEvaluatorFactory<ScalarT>::buildSTKMeshFactory(const Teuchos::ParameterList & mesh_params) const
  {
    Teuchos::RCP<panzer_stk::STK_MeshFactory> mesh_factory;

    // first contruct the mesh factory
    if (mesh_params.get<std::string>("Method") ==  "Exodus File") {
      mesh_factory = Teuchos::rcp(new panzer_stk::STK_ExodusReaderFactory());
      mesh_factory->setParameterList(Teuchos::rcp(new Teuchos::ParameterList(mesh_params.sublist("Exodus File"))));
    }
    else if (mesh_params.get<std::string>("Method") ==  "Pamgen Mesh") {
      mesh_factory = Teuchos::rcp(new panzer_stk::STK_ExodusReaderFactory());
      Teuchos::RCP<Teuchos::ParameterList> pamgenList = Teuchos::rcp(new Teuchos::ParameterList(mesh_params.sublist("Pamgen Mesh")));
      pamgenList->set("File Type","Pamgen"); // For backwards compatibility when pamgen had separate factory from exodus
      mesh_factory->setParameterList(pamgenList);
    }
    else if (mesh_params.get<std::string>("Method") ==  "CGNS File") {
      mesh_factory = Teuchos::rcp(new panzer_stk::STK_ExodusReaderFactory());
      Teuchos::RCP<Teuchos::ParameterList> cgnsList = Teuchos::rcp(new Teuchos::ParameterList(mesh_params.sublist("CGNS File")));
      cgnsList->set("File Type","cgns"); // For backwards compatibility when pamgen had separate factory from exodus
      mesh_factory->setParameterList(cgnsList);
    }
    else if (mesh_params.get<std::string>("Method") ==  "Inline Mesh") {

      int dimension = mesh_params.sublist("Inline Mesh").get<int>("Mesh Dimension");
      std::string typeStr = "";
      if(mesh_params.sublist("Inline Mesh").isParameter("Type"))
         typeStr = mesh_params.sublist("Inline Mesh").get<std::string>("Type");

      if (dimension == 1) {
        mesh_factory = Teuchos::rcp(new panzer_stk::LineMeshFactory);
        Teuchos::RCP<Teuchos::ParameterList> in_mesh = Teuchos::rcp(new Teuchos::ParameterList);
        *in_mesh = mesh_params.sublist("Inline Mesh").sublist("Mesh Factory Parameter List");
        mesh_factory->setParameterList(in_mesh);
      }
      else if (dimension == 2 && typeStr=="Tri") {
        mesh_factory = Teuchos::rcp(new panzer_stk::SquareTriMeshFactory);
        Teuchos::RCP<Teuchos::ParameterList> in_mesh = Teuchos::rcp(new Teuchos::ParameterList);
        *in_mesh = mesh_params.sublist("Inline Mesh").sublist("Mesh Factory Parameter List");
        mesh_factory->setParameterList(in_mesh);
      }
      else if (dimension == 2) {
        mesh_factory = Teuchos::rcp(new panzer_stk::SquareQuadMeshFactory);
        Teuchos::RCP<Teuchos::ParameterList> in_mesh = Teuchos::rcp(new Teuchos::ParameterList);
        *in_mesh = mesh_params.sublist("Inline Mesh").sublist("Mesh Factory Parameter List");
        mesh_factory->setParameterList(in_mesh);
      }
      else if (dimension == 3 && typeStr=="Tet") {
        mesh_factory = Teuchos::rcp(new panzer_stk::CubeTetMeshFactory);
        Teuchos::RCP<Teuchos::ParameterList> in_mesh = Teuchos::rcp(new Teuchos::ParameterList);
        *in_mesh = mesh_params.sublist("Inline Mesh").sublist("Mesh Factory Parameter List");
        mesh_factory->setParameterList(in_mesh);
      }
      else if(dimension == 3) {
        mesh_factory = Teuchos::rcp(new panzer_stk::CubeHexMeshFactory);
        Teuchos::RCP<Teuchos::ParameterList> in_mesh = Teuchos::rcp(new Teuchos::ParameterList);
        *in_mesh = mesh_params.sublist("Inline Mesh").sublist("Mesh Factory Parameter List");
        mesh_factory->setParameterList(in_mesh);
      }
      else if(dimension==4) { // not really "dimension==4" simply a flag to try this other mesh for testing
        mesh_factory = Teuchos::rcp(new panzer_stk::MultiBlockMeshFactory);
        Teuchos::RCP<Teuchos::ParameterList> in_mesh = Teuchos::rcp(new Teuchos::ParameterList);
        *in_mesh = mesh_params.sublist("Inline Mesh").sublist("Mesh Factory Parameter List");
        mesh_factory->setParameterList(in_mesh);
      }
    }
    else if (mesh_params.get<std::string>("Method") ==  "Custom Mesh") {
      mesh_factory = Teuchos::rcp(new panzer_stk::CustomMeshFactory());
      mesh_factory->setParameterList(Teuchos::rcp(new Teuchos::ParameterList(mesh_params.sublist("Custom Mesh"))));
    }
    else {
      // throw a runtime exception for invalid parameter values
    }


    // get rebalancing parameters
    if(mesh_params.isSublist("Rebalance")) {
      const Teuchos::ParameterList & rebalance = mesh_params.sublist("Rebalance");

      // check to see if its enabled
      bool enabled = false;
      if(rebalance.isType<bool>("Enabled"))
        enabled = rebalance.get<bool>("Enabled");

      // we can also use a list description of what to load balance
      Teuchos::RCP<Teuchos::ParameterList> rebalanceCycles;
      if(enabled && rebalance.isSublist("Cycles"))
        rebalanceCycles = Teuchos::rcp(new Teuchos::ParameterList(rebalance.sublist("Cycles")));

      // setup rebalancing as neccessary
      mesh_factory->enableRebalance(enabled,rebalanceCycles);
    }

    return mesh_factory;
  }

  template<typename ScalarT>
  void ModelEvaluatorFactory<ScalarT>::finalizeMeshConstruction(const panzer_stk::STK_MeshFactory & mesh_factory,
                                                                const std::vector<Teuchos::RCP<panzer::PhysicsBlock> > & physicsBlocks,
                                                                const Teuchos::MpiComm<int> mpi_comm,
                                                                panzer_stk::STK_Interface & mesh) const
  {
    // finish building mesh, set required field variables and mesh bulk data
    {
      std::vector<Teuchos::RCP<panzer::PhysicsBlock> >::const_iterator physIter;
      for(physIter=physicsBlocks.begin();physIter!=physicsBlocks.end();++physIter) {
        // what is the block weight for this element block?
        double blockWeight = 0.0;

        Teuchos::RCP<const panzer::PhysicsBlock> pb = *physIter;
        const std::vector<panzer::StrPureBasisPair> & blockFields = pb->getProvidedDOFs();
        const std::vector<std::vector<std::string> > & coordinateDOFs = pb->getCoordinateDOFs();
          // these are treated specially

        // insert all fields into a set
        std::set<panzer::StrPureBasisPair,panzer::StrPureBasisComp> fieldNames;
        fieldNames.insert(blockFields.begin(),blockFields.end());

        // Now we will set up the coordinate fields (make sure to remove
        // the DOF fields)
        {
          std::set<std::string> fields_to_remove;

          // add mesh coordinate fields, setup their removal from fieldNames
          // set to prevent duplication
          for(std::size_t i=0;i<coordinateDOFs.size();i++) {
            mesh.addMeshCoordFields(pb->elementBlockID(),coordinateDOFs[i],"DISPL");
            for(std::size_t j=0;j<coordinateDOFs[i].size();j++)
              fields_to_remove.insert(coordinateDOFs[i][j]);
          }

          // remove the already added coordinate fields
          std::set<std::string>::const_iterator rmItr;
          for (rmItr=fields_to_remove.begin();rmItr!=fields_to_remove.end();++rmItr)
            fieldNames.erase(fieldNames.find(panzer::StrPureBasisPair(*rmItr,Teuchos::null)));
        }

        // add basis to DOF manager: block specific
        std::set<panzer::StrPureBasisPair,panzer::StrPureBasisComp>::const_iterator fieldItr;
        for (fieldItr=fieldNames.begin();fieldItr!=fieldNames.end();++fieldItr) {

          if(fieldItr->second->isScalarBasis() &&
             fieldItr->second->getElementSpace()==panzer::PureBasis::CONST) {
             mesh.addCellField(fieldItr->first,pb->elementBlockID());
          }
          else if(fieldItr->second->isScalarBasis()) {
             mesh.addSolutionField(fieldItr->first,pb->elementBlockID());
          }
          else if(fieldItr->second->isVectorBasis()) {
            std::string d_mod[3] = {"X","Y","Z"};
            for(int d=0;d<fieldItr->second->dimension();d++)
              mesh.addCellField(fieldItr->first+d_mod[d],pb->elementBlockID());
          }
          else { TEUCHOS_ASSERT(false); }

          blockWeight += double(fieldItr->second->cardinality());
        }

        // set the compute block weight (this is the sum of the cardinality of all basis
        // functions on this block
        mesh.setBlockWeight(pb->elementBlockID(),blockWeight);
      }

      mesh_factory.completeMeshConstruction(mesh,*(mpi_comm.getRawMpiComm()));
    }
  }


  template<typename ScalarT>
  Teuchos::RCP<Thyra::ModelEvaluator<ScalarT> > ModelEvaluatorFactory<ScalarT>::getPhysicsModelEvaluator()
  {
    TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(m_physics_me), std::runtime_error,
                       "Objects are not built yet!  Please call buildObjects() member function.");
    return  m_physics_me;
  }

  template<typename ScalarT>
  void ModelEvaluatorFactory<ScalarT>::setUserWorksetFactory(Teuchos::RCP<panzer_stk::WorksetFactory>& user_wkst_factory)
  {
    m_user_wkst_factory = user_wkst_factory;
  }

  template<typename ScalarT>
  Teuchos::RCP<Thyra::ModelEvaluator<ScalarT> > ModelEvaluatorFactory<ScalarT>::getResponseOnlyModelEvaluator()
  {
    if(m_rome_me==Teuchos::null)
      m_rome_me = buildControlModelEvaluator(m_physics_me,m_global_data);

    return m_rome_me;
  }

  template<typename ScalarT>
  Teuchos::RCP<Thyra::ModelEvaluator<ScalarT> > ModelEvaluatorFactory<ScalarT>::
  buildControlModelEvaluator(const Teuchos::RCP<Thyra::ModelEvaluator<ScalarT> > & thyra_me,
                                   const Teuchos::RCP<panzer::GlobalData>& global_data )
  {
    TEUCHOS_TEST_FOR_EXCEPTION(is_null(m_lin_obj_factory), std::runtime_error,
                       "Objects are not built yet!  Please call buildObjects() member function.");
    TEUCHOS_TEST_FOR_EXCEPTION(is_null(m_global_indexer), std::runtime_error,
                       "Objects are not built yet!  Please call buildObjects() member function.");
    TEUCHOS_TEST_FOR_EXCEPTION(is_null(m_mesh), std::runtime_error,
                       "Objects are not built yet!  Please call buildObjects() member function.");

    bool useCoordinateUpdate = false;

    Teuchos::ParameterList& p = *this->getNonconstParameterList();
    
    Teuchos::RCP<NOXObserverFactory> nox_observer_factory = 
        Teuchos::rcp(new NOXObserverFactory(m_stkIOResponseLibrary));
    Teuchos::RCP<Teuchos::ParameterList> observers_to_build = 
	          Teuchos::parameterList(p.sublist("Visual Output"));
	  nox_observer_factory->setParameterList(observers_to_build);

    Teuchos::ParameterList & solncntl_params = p.sublist("Solution Control");
    Teuchos::RCP<Teuchos::ParameterList> piro_params = Teuchos::rcp(new Teuchos::ParameterList(solncntl_params));
    Teuchos::RCP<Thyra::ModelEvaluatorDefaultBase<double> > piro;

    std::string solver = solncntl_params.get<std::string>("Piro Solver");
    Teuchos::RCP<Thyra::ModelEvaluatorDefaultBase<double> > thyra_me_db
       = Teuchos::rcp_dynamic_cast<Thyra::ModelEvaluatorDefaultBase<double> >(thyra_me);

    if ( solver=="NOX" ) {
      Teuchos::RCP<NOX::Abstract::PrePostOperator> ppo = nox_observer_factory->buildNOXObserver(m_mesh,m_global_indexer,m_lin_obj_factory);
      piro_params->sublist("NOX").sublist("Solver Options").set("User Defined Pre/Post Operator", ppo);

      piro = Teuchos::rcp(new Piro::NOXSolver<double>(piro_params,thyra_me_db));
      TEUCHOS_ASSERT(nonnull(piro));

      // override printing to use panzer ostream
      piro_params->sublist("NOX").sublist("Printing").set<Teuchos::RCP<std::ostream> >("Output Stream",global_data->os);
      piro_params->sublist("NOX").sublist("Printing").set<Teuchos::RCP<std::ostream> >("Error Stream",global_data->os);
      piro_params->sublist("NOX").sublist("Printing").set<int>("Output Processor",global_data->os->getOutputToRootOnly());
    }

    else if ( solver == "LOCA" ) {
      const Teuchos::RCP<LOCA::Thyra::SaveDataStrategy> saveEigData =
        Teuchos::rcp( new LOCA::Thyra::SaveDataStrategy() );

      Teuchos::RCP<NOX::Abstract::PrePostOperator> ppo = nox_observer_factory->buildNOXObserver(m_mesh,m_global_indexer,m_lin_obj_factory);
      piro_params->sublist("NOX").sublist("Solver Options").set("User Defined Pre/Post Operator", ppo);

      piro = Teuchos::rcp(new Piro::LOCASolver<double>(piro_params,thyra_me_db,Teuchos::null,saveEigData));
      TEUCHOS_ASSERT(nonnull(piro));

      // override printing to use panzer ostream
      piro_params->sublist("NOX").sublist("Printing").set<Teuchos::RCP<std::ostream> >("Output Stream",global_data->os);
      piro_params->sublist("NOX").sublist("Printing").set<Teuchos::RCP<std::ostream> >("Error Stream",global_data->os);
      piro_params->sublist("NOX").sublist("Printing").set<int>("Output Processor",global_data->os->getOutputToRootOnly());
    }
    else if (solver=="Rythmos") {

      Teuchos::RCP<const RythmosObserverFactory> observer_factory =
        Teuchos::rcp(new RythmosObserverFactory(m_stkIOResponseLibrary,m_response_library->getWorksetContainer(),useCoordinateUpdate));

      // install the nox observer
      if(observer_factory->useNOXObserver()) {
        Teuchos::RCP<NOX::Abstract::PrePostOperator> ppo = nox_observer_factory->buildNOXObserver(m_mesh,m_global_indexer,m_lin_obj_factory);
        piro_params->sublist("NOX").sublist("Solver Options").set("User Defined Pre/Post Operator", ppo);
      }

      // override printing to use panzer ostream
      piro_params->sublist("NOX").sublist("Printing").set<Teuchos::RCP<std::ostream> >("Output Stream",global_data->os);
      piro_params->sublist("NOX").sublist("Printing").set<Teuchos::RCP<std::ostream> >("Error Stream",global_data->os);
      piro_params->sublist("NOX").sublist("Printing").set<int>("Output Processor",global_data->os->getOutputToRootOnly());

      // use the user specfied rythmos solver if they pass one in
      Teuchos::RCP<Piro::RythmosSolver<double> > piro_rythmos = Teuchos::rcp(new Piro::RythmosSolver<double>());

      // if you are using explicit RK, make sure to wrap the ME in an explicit model evaluator decorator
      Teuchos::RCP<Thyra::ModelEvaluator<ScalarT> > rythmos_me = thyra_me;
      const std::string stepper_type = piro_params->sublist("Rythmos").get<std::string>("Stepper Type");
      if(stepper_type=="Explicit RK" || stepper_type=="Forward Euler") {
        const Teuchos::ParameterList & assembly_params = p.sublist("Assembly");
        bool lumpExplicitMass = assembly_params.get<bool>("Lump Explicit Mass");
        rythmos_me = Teuchos::rcp(new panzer::ExplicitModelEvaluator<ScalarT>(thyra_me,!useDynamicCoordinates_,lumpExplicitMass));
      }

      piro_rythmos->initialize(piro_params, rythmos_me, observer_factory->buildRythmosObserver(m_mesh,m_global_indexer,m_lin_obj_factory));

      piro = piro_rythmos;
    }
    else if (solver=="Tempus") {
    //  Teuchos::RCP<Thyra::ModelEvaluator<ScalarT> > tempus_me = thyra_me;
    //  piro = Teuchos::rcp(new Piro::TempusSolver<double>(piro_params,thyra_me));
      const Teuchos::RCP<Teuchos::ParameterList> & pl = sublist(piro_params, "Tempus");
  //    pl->print(std::cout);
      Piro::SENS_METHOD sens_method; 
      std::string sens_method_string = piro_params->get<std::string>("Sensitivity Method", "None");
      if (sens_method_string == "None") sens_method = Piro::NONE; 
      else if (sens_method_string == "Forward") sens_method = Piro::FORWARD; 
      else if (sens_method_string == "Adjoint") sens_method = Piro::ADJOINT; 
      Teuchos::RCP<Piro::TempusIntegrator<double> > integrator 
          = Teuchos::rcp(new Piro::TempusIntegrator<double>(pl, thyra_me_db, sens_method));

      // observer
      Teuchos::RCP<const TempusObserverFactory> observer_factory =
        Teuchos::rcp(new TempusObserverFactory(m_stkIOResponseLibrary,m_response_library->getWorksetContainer(),useCoordinateUpdate));
      int freq = observers_to_build->get<int>("Output Frequency");
      Teuchos::RCP< Tempus::IntegratorObserver<double> > observer = 
        observer_factory->buildTempusObserver(m_mesh,m_global_indexer,m_lin_obj_factory,freq);
      integrator -> setObserver(observer);

      Teuchos::RCP<Thyra::NonlinearSolverBase<double> > stepSolver;
      if( piro_params->isSublist("NOX")) {
        stepSolver =  Teuchos::rcp(new Thyra::NOXNonlinearSolver);
        Teuchos::RCP<Teuchos::ParameterList> nox_params = sublist(piro_params, "NOX");
        stepSolver->setParameterList(nox_params);
      }
      else
        stepSolver = Teuchos::null;

      Teuchos::RCP<Teuchos::ParameterList> stepperPL = Teuchos::rcp(&(pl->sublist("Tempus Stepper")), false);

      Teuchos::RCP<Tempus::StepperFactory<double> > sf = Teuchos::rcp(new Tempus::StepperFactory<double>());
      const Teuchos::RCP<Tempus::Stepper<double> > stepper = sf->createStepper(stepperPL, thyra_me_db);
      if( stepper->isExplicit() )
      {}
    /*  Teuchos::RCP<const Thyra::ModelEvaluator<double>>  thyra_me_const 
        = Teuchos::rcp_const_cast<const Thyra::ModelEvaluator<double> >( thyra_me );
      auto solutionHistory = Tempus::createSolutionHistoryME(thyra_me_const);
      stepper->setInitialConditions(solutionHistory);*/
      double initialTime = pl->sublist("Tempus Integrator").sublist("Time Step Control").get<double>("Initial Time",0.0);
      double finalTime = pl->sublist("Tempus Integrator").sublist("Time Step Control").get<double>("Final Time");

    /*  auto x0 =  thyra_me_db->getNominalValues().get_x(); 
      std::cout << *x0 << std::endl;
      const int num_param = thyra_me_db->get_p_space(0)->dim();
      Teuchos::RCP<Thyra::MultiVectorBase<double> > DxDp0 =
        Thyra::createMembers(thyra_me_db->get_x_space(), num_param);
      DxDp0->assign(0.0);
      integrator->initializeSolutionHistory(initialTime, x0);*/

      Teuchos::RCP<Piro::TempusSolver<double> > piro_tempus =
        Teuchos::rcp(new Piro::TempusSolver<double>(integrator, stepper, stepSolver, thyra_me_db, finalTime, sens_method_string));
    //  piro_tempus -> setInitialState(initialTime, Teuchos::rcp_const_cast<Thyra::VectorBase< double > >(x0));

      piro = piro_tempus;
    }
    else if (solver=="Eigen") 
    {
      const Teuchos::RCP<Teuchos::ParameterList> pl = sublist(piro_params, "Eigen");
      std::shared_ptr< EigenObserver_WriteToExodus<double> >  sptr = 
        std::make_shared< EigenObserver_WriteToExodus<double> >(m_mesh,m_global_indexer,m_lin_obj_factory,m_stkIOResponseLibrary);
      piro = buildEigenSolver<double>(pl, thyra_me_db, sptr);
      TEUCHOS_ASSERT(nonnull(piro));
    }
	else if (solver=="ROL") {
		Teuchos::RCP<Teuchos::ParameterList> pl = sublist(piro_params, "ROL");
		piro = Teuchos::rcp(new ROLSolver<double>(*pl, thyra_me_db));
	}
    else {
      TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
                         "Error: Unknown Piro Solver : " << solver);
    }
    return piro;
  }

  template<typename ScalarT>
  Teuchos::RCP<FieldManagerBuilder>
  ModelEvaluatorFactory<ScalarT>::
  buildFieldManagerBuilder(const Teuchos::RCP<panzer::WorksetContainer> & wc,
                           const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
                           const std::vector<panzer::BC> & bcs,
                           const std::vector< Dirichlet<panzer::Traits::RealType> >& ds,
                           const std::vector< NeumannBoundary >& ns,
                           const std::vector< CLoad >& cloads,
                           const panzer::EquationSetFactory & eqset_factory,
                           const panzer::BCStrategyFactory& bc_factory,
                           const panzer::ClosureModelFactory_TemplateManager<panzer::Traits>& volume_cm_factory,
                           const panzer::ClosureModelFactory_TemplateManager<panzer::Traits>& bc_cm_factory,
                           const Teuchos::ParameterList& closure_models,
                           const panzer::LinearObjFactory<panzer::Traits> & lo_factory,
                           const Teuchos::ParameterList& user_data,
                           bool writeGraph,const std::string & graphPrefix,
			   bool write_field_managers,const std::string & field_manager_prefix) const
  {
    Teuchos::RCP<FieldManagerBuilder> fmb = Teuchos::rcp(new FieldManagerBuilder);
    fmb->setWorksetContainer(wc);
    fmb->buildMaterials(user_data);
    fmb->setupVolumeFieldManagers(m_mesh, physicsBlocks,volume_cm_factory,closure_models,lo_factory,user_data);
    fmb->setupBCFieldManagers(bcs,physicsBlocks,eqset_factory,bc_cm_factory,bc_factory,closure_models,lo_factory,user_data);
    fmb->setupDirichletFieldManager(ds, m_mesh);
    fmb->setupNeumannFieldManager(ns,cloads,lo_factory,physicsBlocks);

    // Print Phalanx DAGs
    if (writeGraph){
      fmb->writeVolumeGraphvizDependencyFiles(graphPrefix, physicsBlocks);
      fmb->writeBCGraphvizDependencyFiles(graphPrefix);
    }
    if (write_field_managers){
      fmb->writeVolumeTextDependencyFiles(graphPrefix, physicsBlocks);
      fmb->writeBCTextDependencyFiles(field_manager_prefix);
    }
    
    return fmb;
  }

  // Setup STK response library for writing out the solution fields
  ////////////////////////////////////////////////////////////////////////
  template<typename ScalarT>
  Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > ModelEvaluatorFactory<ScalarT>::
  initializeSolnWriterResponseLibrary(const Teuchos::RCP<panzer::WorksetContainer> & wc,
                                      const Teuchos::RCP<const panzer::GlobalIndexer> & ugi,
                                      const Teuchos::RCP<const panzer::LinearObjFactory<panzer::Traits> > & lof,
                                      const Teuchos::RCP<panzer_stk::STK_Interface> & mesh) const
  {
     Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > stkIOResponseLibrary
        = Teuchos::rcp(new panzer::ResponseLibrary<panzer::Traits>(wc,ugi,lof));

     std::vector<std::string> eBlocks;
     mesh->getElementBlockNames(eBlocks);

     panzer_stk::RespFactorySolnWriter_Builder builder;
     builder.mesh = mesh;
     stkIOResponseLibrary->addResponse("Main Field Output",eBlocks,builder);

     return stkIOResponseLibrary;
  }

  template<typename ScalarT>
  void ModelEvaluatorFactory<ScalarT>::
  finalizeSolnWriterResponseLibrary(panzer::ResponseLibrary<panzer::Traits> & rl,
                                    const std::vector<Teuchos::RCP<panzer::PhysicsBlock> > & physicsBlocks,
                                    const panzer::ClosureModelFactory_TemplateManager<panzer::Traits> & cm_factory,
                                    const Teuchos::ParameterList & closure_models,
                                    int workset_size, Teuchos::ParameterList & user_data) const
  {
     user_data.set<int>("Workset Size",workset_size);
     rl.buildResponseEvaluators(physicsBlocks, cm_factory, closure_models, user_data);
  }

  template<typename ScalarT>
  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> > ModelEvaluatorFactory<ScalarT>::
  buildLOWSFactory(bool blockedAssembly,
                   const Teuchos::RCP<const panzer::GlobalIndexer> & globalIndexer,
                   const Teuchos::RCP<panzer::ConnManager> & conn_manager,
                   const Teuchos::RCP<panzer_stk::STK_Interface> & mesh,
                   const Teuchos::RCP<const Teuchos::MpiComm<int> > & mpi_comm,
                   const Teuchos::RCP<Teko::RequestHandler> & reqHandler
                   ) const
  {
    const Teuchos::ParameterList & p = *this->getParameterList();
    Teuchos::RCP<Teuchos::ParameterList> strat_params = Teuchos::rcp(new Teuchos::ParameterList(p.sublist("Linear Solver")) );

    bool writeCoordinates = false;
    if(p.sublist("Debug Options").isType<bool>("Write Coordinates"))
      writeCoordinates = p.sublist("Debug Options").get<bool>("Write Coordinates");

    bool writeTopo = false;
    if(p.sublist("Debug Options").isType<bool>("Write Topology"))
      writeTopo = p.sublist("Debug Options").get<bool>("Write Topology");


    return panzer_stk::buildLOWSFactory(
                            blockedAssembly,globalIndexer,conn_manager,
                            Teuchos::as<int>(mesh->getDimension()), mpi_comm, strat_params,
                            reqHandler,
                            writeCoordinates,
                            writeTopo );
  }

  template<typename ScalarT>
  void ModelEvaluatorFactory<ScalarT>::
  buildResponses(const panzer::ClosureModelFactory_TemplateManager<panzer::Traits> & cm_factory,
                 const bool write_graphviz_file,
                 const std::string& graphviz_file_prefix)
  {
    typedef panzer::ModelEvaluator<double> PanzerME;
    Teuchos::RCP<PanzerME> panzer_me = Teuchos::rcp_dynamic_cast<PanzerME>(m_physics_me);

    Teuchos::RCP<const Teuchos::ParameterList> input_params = getMyParamList();
    Teuchos::ParameterList responses = input_params->sublist("Responses");

    // add a volume response functional for each field 
    for(Teuchos::ParameterList::ConstIterator itr=responses.begin();itr!=responses.end();++itr) {
        const std::string name = responses.name(itr);
        TEUCHOS_ASSERT(responses.entry(itr).isList());
        Teuchos::ParameterList & lst = Teuchos::getValue<Teuchos::ParameterList>(responses.entry(itr));

        // parameterize the builder
        panzer::FunctionalResponse_Builder<int,int> builder;
        builder.comm = MPI_COMM_WORLD; // good enough
        builder.cubatureDegree = 2;
        builder.requiresCellIntegral = lst.isType<bool>("Requires Cell Integral") ? lst.get<bool>("Requires Cell Integral"): false;
        builder.quadPointField = lst.get<std::string>("Field Name");

        // add the respone
        std::vector<std::string> eblocks;
        panzer::StringTokenizer(eblocks,lst.get<std::string>("Element Blocks"),",",true);
        
        std::vector<panzer::WorksetDescriptor> wkst_descs;
        for(std::size_t i=0;i<eblocks.size();i++) 
          wkst_descs.push_back(panzer::blockDescriptor(eblocks[i]));

        int respIndex = panzer_me->addResponse(name,wkst_descs,builder);
      //  panzer_me->addResponse(responseName,wkstDesc,builder);
        m_responseIndexToName[respIndex] = name;
    }
 

    Teuchos::ParameterList & p = *this->getNonconstParameterList();
    Teuchos::ParameterList & user_data = p.sublist("User Data");
    Teuchos::ParameterList & closure_models = p.sublist("Closure Models");

    if(panzer_me!=Teuchos::null) {
      panzer_me->buildResponses(m_physics_blocks,*m_eqset_factory,cm_factory,closure_models,user_data,write_graphviz_file,graphviz_file_prefix);
      return;
    }

    TEUCHOS_ASSERT(false);
  }

  template<typename ScalarT>
  void ModelEvaluatorFactory<ScalarT>::
  buildIOResponses(const panzer::ClosureModelFactory_TemplateManager<panzer::Traits> & cm_factory,
                  Teuchos::RCP<Teuchos::ParameterList>& input_params,
                 const bool write_graphviz_file,
                 const std::string& graphviz_file_prefix)
  {
  //  m_stkIOResponseLibrary = Teuchos::rcp(new panzer::ResponseLibrary<panzer::Traits>());
    m_stkIOResponseLibrary->initialize(*m_response_library);

    Teuchos::ParameterList user_data(input_params->sublist("User Data"));
    user_data.set<int>("Workset Size",input_params->sublist("Assembly").get<int>("Workset Size"));
    
    m_stkIOResponseLibrary->buildResponseEvaluators(m_physics_blocks,
                                        cm_factory,
                                        input_params->sublist("Closure Models"),
                                        user_data);
  }

  template<typename ScalarT>
  void ModelEvaluatorFactory<ScalarT>::
  buildFluxResponses(const panzer::ClosureModelFactory_TemplateManager<panzer::Traits> & cm_factory,
                  Teuchos::RCP<Teuchos::ParameterList>& input_params,
                 const bool write_graphviz_file,
                 const std::string& graphviz_file_prefix)
  {
    m_fluxResponse = Teuchos::rcp(new panzer::ResponseLibrary<panzer::Traits>());
    m_fluxResponse->initialize(*m_response_library);

    if( !input_params->isSublist("fluxResponses") ) return;
    Teuchos::ParameterList responses = input_params->sublist("fluxResponses");
    if( responses.numParams()<=0 ) return;

    // build high-order flux response
    for(Teuchos::ParameterList::ConstIterator itr=responses.begin();itr!=responses.end();++itr) {
        const std::string name = responses.name(itr);
        TEUCHOS_ASSERT(responses.entry(itr).isList());
        Teuchos::ParameterList & lst = Teuchos::getValue<Teuchos::ParameterList>(responses.entry(itr));

        SiYuan::Response_Builder builder;
        builder.comm = MPI_COMM_WORLD;
        builder.cubatureDegree = 2;
        builder.dof_name = lst.get<std::string>("DOF name");
        TEUCHOS_ASSERT( builder.dof_name.length()>0 );
  
        std::vector<panzer::WorksetDescriptor> sidesets;
        std::string eblock = lst.get<std::string>("Element Block");
        std::string side = lst.get<std::string>("Sideset ID");
        sidesets.push_back(panzer::sidesetVolumeDescriptor(eblock,side));
  
        m_fluxResponse->addResponse(name,sidesets,builder);
    }
  
      {
        Teuchos::ParameterList user_data(input_params->sublist("User Data"));
        user_data.set<int>("Workset Size",input_params->sublist("Assembly").get<int>("Workset Size"));
      
        m_fluxResponse->buildResponseEvaluators(m_physics_blocks,
                                          *m_eqset_factory,
                                          cm_factory,
                                          input_params->sublist("Closure Models"),
                                          user_data);
      }

      for(Teuchos::ParameterList::ConstIterator itr=responses.begin();itr!=responses.end();++itr) {
        const std::string name = responses.name(itr);
        Teuchos::RCP<panzer::ResponseMESupportBase<panzer::Traits::Residual> > resp
            = Teuchos::rcp_dynamic_cast<panzer::ResponseMESupportBase<panzer::Traits::Residual> >(m_fluxResponse->getResponse<panzer::Traits::Residual>(name),true);

        const Teuchos::RCP<Thyra::VectorBase<double> >  vec = Thyra::createMember(resp->getVectorSpace());
        resp->setVector( vec );
      }
  }

  template<typename ScalarT>
  void ModelEvaluatorFactory<ScalarT>::
  evalResponses(Teuchos::RCP<Thyra::VectorBase<double> >& gx)
  {
    Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::rcp(new Teuchos::FancyOStream(Teuchos::rcp(&std::cout,false)));
    if(m_physics_me->Ng()>0) {

        Thyra::ModelEvaluatorBase::InArgs<double> respInArgs = m_physics_me->createInArgs();
        Thyra::ModelEvaluatorBase::OutArgs<double> respOutArgs = m_physics_me->createOutArgs();

        TEUCHOS_ASSERT(m_physics_me->Ng()==respOutArgs.Ng());
   
        respInArgs.set_x(gx);
   
        // set up response out args
        for(int i=0;i<respOutArgs.Ng();i++) {
	        Teuchos::RCP<Thyra::VectorBase<double> > response = Thyra::createMember(*m_physics_me->get_g_space(i));
          respOutArgs.set_g(i,response);
        }
   
        // Now, solve the problem and return the responses
        m_physics_me->evalModel(respInArgs, respOutArgs);
   
        // loop over out args for printing
        for(int i=0;i<respOutArgs.Ng();i++) {
	        Teuchos::RCP<Thyra::VectorBase<double> > response = respOutArgs.get_g(i);

          TEUCHOS_ASSERT(response!=Teuchos::null); // should not be null!

            *out << "Response Value \"" << m_responseIndexToName[i] << "\": " << *response << std::endl;
        }
    }
    if( m_fluxResponse->responseEvaluatorsBuilt() ) {
      panzer::AssemblyEngineInArgs ae_inargs;
      ae_inargs.container_ = m_lin_obj_factory->buildLinearObjContainer();
      ae_inargs.ghostedContainer_ = m_lin_obj_factory->buildGhostedLinearObjContainer();
      ae_inargs.alpha = 0.0;
      ae_inargs.beta = 1.0;
      ae_inargs.evaluate_transient_terms = false;

      // initialize the ghosted container
      m_lin_obj_factory->initializeGhostedContainer(panzer::LinearObjContainer::X,*ae_inargs.ghostedContainer_);

      if( m_blockedAssembly ) {
        const Teuchos::RCP<panzer::BlockedTpetraLinearObjContainer<double,int,panzer::GlobalOrdinal>> tpGlobalContainer
           = Teuchos::rcp_dynamic_cast<panzer::BlockedTpetraLinearObjContainer<double,int,panzer::GlobalOrdinal>>(ae_inargs.container_,true);
        tpGlobalContainer->set_x_th(gx);
      } else {
        const Teuchos::RCP<panzer::TpetraLinearObjContainer<double,int,panzer::GlobalOrdinal>> tpGlobalContainer
           = Teuchos::rcp_dynamic_cast<panzer::TpetraLinearObjContainer<double,int,panzer::GlobalOrdinal>>(ae_inargs.container_,true);
        tpGlobalContainer->set_x_th(gx);
      }

      // evaluate current on contacts
      m_fluxResponse->addResponsesToInArgs<panzer::Traits::Residual>(ae_inargs);
      m_fluxResponse->evaluate<panzer::Traits::Residual>(ae_inargs);

      std::vector<Teuchos::RCP<panzer::ResponseBase> > v;
      m_fluxResponse->getResponses<panzer::Traits::Residual>(v);

        // output current values
      *out << "\nFlux values: \n";
      for(auto r: v ) {
        std::string currentRespName = r->getName();

        Teuchos::RCP<panzer::Response_Functional<panzer::Traits::Residual> > resp
            = Teuchos::rcp_dynamic_cast<panzer::Response_Functional<panzer::Traits::Residual> >(m_fluxResponse->getResponse<panzer::Traits::Residual>(currentRespName),true);

        *out << "   " << currentRespName << " = " << resp->value << std::endl;
      }
    }
  }
}

#endif
